#Code for "Evaluating Punctuated Equilibrium Dynamics Within a Crisis Context"
#Saahir Shafi & Daniel J. Mallinson
#07-03-2023


rm(list=ls())


#install.packages(c("foreign", "moments", "Hmisc", "EnvStats", "lmom"))

library(foreign)
library(moments)
library(Hmisc)
library(EnvStats)
library(lmom)

## Calculate Kurtosis by System and Country_Weekly

weekly <- read.csv("PET_weekly_5.2023.csv")

weekly <- na.omit(weekly)

#Country

countries <- unique(weekly[,1])

for(i in 1:length(countries)){
  use.data <- weekly[which(weekly[,1]==countries[i]),]
  #Stringency
  l.moment <- round(samlmu(use.data$Percent_changes_in_SI), digits=2)
  l.kurtosis.si <- l.moment[4]
  #GRI
  l.moment <- round(samlmu(use.data$Percent_changes_in_GRI), digits=2)
  l.kurtosis.gri <- l.moment[4]
  #CHI
  l.moment <- round(samlmu(use.data$Percent_changes_in_CHI), digits=2)
  l.kurtosis.chi <- l.moment[4]
  #ESI
  l.moment <- round(samlmu(use.data$Percent_changes_in_ESI), digits=2)
  l.kurtosis.esi <- l.moment[4]
  if(i==1){
    country.kurt <- as.data.frame(cbind(countries[i], l.kurtosis.si, l.kurtosis.gri, l.kurtosis.chi, l.kurtosis.esi, unique(use.data$electoral), unique(use.data$libdem), unique(use.data$partipdem), unique(use.data$delibdem), unique(use.data$freedomhouse)))
  }else{
    bind <- cbind(countries[i], l.kurtosis.si, l.kurtosis.gri, l.kurtosis.chi, l.kurtosis.esi, unique(use.data$electoral), unique(use.data$libdem), unique(use.data$partipdem), unique(use.data$delibdem), unique(use.data$freedomhouse))
    country.kurt <- rbind(country.kurt, bind)
  }
}

write.csv(country.kurt, file="kurtosis_by_country_weekly.csv")

#Government System - Electoral Democracy
system.electoral <- unique(weekly$electoral)

for(i in 1:length(system.electoral)){
  use.data <- weekly[which(weekly$electoral==system.electoral[i]),]
  #Stringency
  l.moment <- round(samlmu(use.data$Percent_changes_in_SI), digits=2)
  l.kurtosis.si <- l.moment[4]
  #GRI
  l.moment <- round(samlmu(use.data$Percent_changes_in_GRI), digits=2)
  l.kurtosis.gri <- l.moment[4]
  #CHI
  l.moment <- round(samlmu(use.data$Percent_changes_in_CHI), digits=2)
  l.kurtosis.chi <- l.moment[4]
  #ESI
  l.moment <- round(samlmu(use.data$Percent_changes_in_ESI), digits=2)
  l.kurtosis.esi <- l.moment[4]
  if(i==1){
    system.kurt.electoral <- as.data.frame(cbind(system.electoral[i], l.kurtosis.si, l.kurtosis.gri, l.kurtosis.chi, l.kurtosis.esi))
  }else{
    bind <- cbind(system.electoral[i], l.kurtosis.si, l.kurtosis.gri, l.kurtosis.chi, l.kurtosis.esi)
    system.kurt.electoral <- rbind(system.kurt.electoral, bind)
  }
}

write.csv(system.kurt.electoral, file="kurtosis_by_system_weekly_electoral.csv")


#Government System - Liberal Democracy
system.libdem <- unique(weekly$libdem)

for(i in 1:length(system.libdem)){
  use.data <- weekly[which(weekly$libdem==system.libdem[i]),]
  #Stringency
  l.moment <- round(samlmu(use.data$Percent_changes_in_SI), digits=2)
  l.kurtosis.si <- l.moment[4]
  #GRI
  l.moment <- round(samlmu(use.data$Percent_changes_in_GRI), digits=2)
  l.kurtosis.gri <- l.moment[4]
  #CHI
  l.moment <- round(samlmu(use.data$Percent_changes_in_CHI), digits=2)
  l.kurtosis.chi <- l.moment[4]
  #ESI
  l.moment <- round(samlmu(use.data$Percent_changes_in_ESI), digits=2)
  l.kurtosis.esi <- l.moment[4]
  if(i==1){
    system.kurt.libdem <- as.data.frame(cbind(system.libdem[i], l.kurtosis.si, l.kurtosis.gri, l.kurtosis.chi, l.kurtosis.esi))
  }else{
    bind <- cbind(system.libdem[i], l.kurtosis.si, l.kurtosis.gri, l.kurtosis.chi, l.kurtosis.esi)
    system.kurt.libdem <- rbind(system.kurt.libdem, bind)
  }
}

write.csv(system.kurt.libdem, file="kurtosis_by_system_weekly_libdem.csv")

#Government System - Participatory Democracy
system.partipdem <- unique(weekly$partipdem)

for(i in 1:length(system.partipdem)){
  use.data <- weekly[which(weekly$partipdem==system.partipdem[i]),]
  #Stringency
  l.moment <- round(samlmu(use.data$Percent_changes_in_SI), digits=2)
  l.kurtosis.si <- l.moment[4]
  #GRI
  l.moment <- round(samlmu(use.data$Percent_changes_in_GRI), digits=2)
  l.kurtosis.gri <- l.moment[4]
  #CHI
  l.moment <- round(samlmu(use.data$Percent_changes_in_CHI), digits=2)
  l.kurtosis.chi <- l.moment[4]
  #ESI
  l.moment <- round(samlmu(use.data$Percent_changes_in_ESI), digits=2)
  l.kurtosis.esi <- l.moment[4]
  if(i==1){
    system.kurt.partipdem <- as.data.frame(cbind(system.partipdem[i], l.kurtosis.si, l.kurtosis.gri, l.kurtosis.chi, l.kurtosis.esi))
  }else{
    bind <- cbind(system.partipdem[i], l.kurtosis.si, l.kurtosis.gri, l.kurtosis.chi, l.kurtosis.esi)
    system.kurt.partipdem <- rbind(system.kurt.partipdem, bind)
  }
}

write.csv(system.kurt.partipdem, file="kurtosis_by_system_weekly_partipdem.csv")

#Government System - Deliberative Democracy
system.delibdem <- unique(weekly$delibdem)

for(i in 1:length(system.delibdem)){
  use.data <- weekly[which(weekly$delibdem==system.delibdem[i]),]
  #Stringency
  l.moment <- round(samlmu(use.data$Percent_changes_in_SI), digits=2)
  l.kurtosis.si <- l.moment[4]
  #GRI
  l.moment <- round(samlmu(use.data$Percent_changes_in_GRI), digits=2)
  l.kurtosis.gri <- l.moment[4]
  #CHI
  l.moment <- round(samlmu(use.data$Percent_changes_in_CHI), digits=2)
  l.kurtosis.chi <- l.moment[4]
  #ESI
  l.moment <- round(samlmu(use.data$Percent_changes_in_ESI), digits=2)
  l.kurtosis.esi <- l.moment[4]
  if(i==1){
    system.kurt.delibdem <- as.data.frame(cbind(system.delibdem[i], l.kurtosis.si, l.kurtosis.gri, l.kurtosis.chi, l.kurtosis.esi))
  }else{
    bind <- cbind(system.delibdem[i], l.kurtosis.si, l.kurtosis.gri, l.kurtosis.chi, l.kurtosis.esi)
    system.kurt.delibdem <- rbind(system.kurt.delibdem, bind)
  }
}

write.csv(system.kurt.delibdem, file="kurtosis_by_system_weekly_delibdem.csv")

#Government System - Freedom House
system.freedomhouse <- unique(weekly$freedomhouse)

for(i in 1:length(system.freedomhouse)){
  use.data <- weekly[which(weekly$freedomhouse==system.freedomhouse[i]),]
  #Stringency
  l.moment <- round(samlmu(use.data$Percent_changes_in_SI), digits=2)
  l.kurtosis.si <- l.moment[4]
  #GRI
  l.moment <- round(samlmu(use.data$Percent_changes_in_GRI), digits=2)
  l.kurtosis.gri <- l.moment[4]
  #CHI
  l.moment <- round(samlmu(use.data$Percent_changes_in_CHI), digits=2)
  l.kurtosis.chi <- l.moment[4]
  #ESI
  l.moment <- round(samlmu(use.data$Percent_changes_in_ESI), digits=2)
  l.kurtosis.esi <- l.moment[4]
  if(i==1){
    system.kurt.freedomhouse <- as.data.frame(cbind(system.freedomhouse[i], l.kurtosis.si, l.kurtosis.gri, l.kurtosis.chi, l.kurtosis.esi))
  }else{
    bind <- cbind(system.freedomhouse[i], l.kurtosis.si, l.kurtosis.gri, l.kurtosis.chi, l.kurtosis.esi)
    system.kurt.freedomhouse <- rbind(system.kurt.freedomhouse, bind)
  }
}

write.csv(system.kurt.freedomhouse, file="kurtosis_by_system_weekly_freedomhouse.csv")


## GRI Kurtosis Distributions

#Load GRI data for kurtosis test
GRI_policy_change <- read.csv("PET_weekly_5.2023.csv")

data <- GRI_policy_change
data <- na.omit(data) #remove blank pct_diff cells

#############################################################################
############################# Figure 4a: Free & GRI Weekly ##################
#############################################################################


##Free

Free.sub <- subset(data, data$freedomhouse_1 == 1)
Free.sub

###
specify_decimal <- function(l,m) format(round(l,m), nsmall=m)

#setEPS()
#postscript("All_Spending_Histogram.eps")
pdf("Policy_Change_Histogram_Free.pdf")
#tiff("Policy_Change_Histogram.tiff", height=6, width=6, units="in", res=600)
x <- Free.sub$Percent_changes_in_GRI
View(x)
kurtosis <- specify_decimal(kurtosis(x), 2)
kurtosis
l.moment <- round(samlmu(x), digits=2)
l.moment
l.kurtosis <- l.moment[4]
h <- hist(x, breaks=50, main="Free Countries", xlab="Percent Changes in GRI (Weekly)", ylab = "Frequency", ylim=c(0,7000), xlim=c(-.6,.8), col="gray50", border="gray50")
xfit <- seq(min(x), max(x), length=100)
yfit <- dnorm(xfit,  mean=mean(x), sd=sd(x))
yfit <- yfit*diff(h$mids[1:2]*length(x))
lines(xfit, yfit, col="black", lwd=2, lty=1)
abline(v=0, col="black", lwd=2)
abline(v=median(x), col="black", lwd=2, lty=2)
kurt.lab <- bquote(bold(L-kurtosis == .(specify_decimal(l.kurtosis,2))))
med.lab <- bquote(bold(Median == .(specify_decimal(median(x),2))))
text(.5, 6000, labels=kurt.lab)
text(.5, 5500, labels=med.lab)
dev.off()

#############################################################################
############################# Figure 4b: Partly Free & GRI Weekly ############
#############################################################################


##Partly Free

PFree.sub <- subset(data, data$freedomhouse_1 == .5)
PFree.sub

###
specify_decimal <- function(l,m) format(round(l,m), nsmall=m)

#setEPS()
#postscript("All_Spending_Histogram.eps")
pdf("Policy_Change_Histogram_PartlyFree.pdf")
#tiff("Policy_Change_Histogram.tiff", height=6, width=6, units="in", res=600)
x <- PFree.sub$Percent_changes_in_GRI
View(x)
kurtosis <- specify_decimal(kurtosis(x), 2)
kurtosis
l.moment <- round(samlmu(x), digits=2)
l.moment
l.kurtosis <- l.moment[4]
h <- hist(x, breaks=50, main="Partly Free Countries", xlab="Percent Changes in GRI (Weekly)", ylab = "Frequency", ylim=c(0,7000), xlim=c(-.6,.8), col="gray50", border="gray50")
xfit <- seq(min(x), max(x), length=100)
yfit <- dnorm(xfit,  mean=mean(x), sd=sd(x))
yfit <- yfit*diff(h$mids[1:2]*length(x))
lines(xfit, yfit, col="black", lwd=2, lty=1)
abline(v=0, col="black", lwd=2)
abline(v=median(x), col="black", lwd=2, lty=2)
kurt.lab <- bquote(bold(L-kurtosis == .(specify_decimal(l.kurtosis,2))))
med.lab <- bquote(bold(Median == .(specify_decimal(median(x),2))))
text(.5, 6000, labels=kurt.lab)
text(.5, 5500, labels=med.lab)
dev.off()

#############################################################################
############################# Figure 4c: Not Free & GRI Weekly ###############
#############################################################################


##Not Free

NFree.sub <- subset(data, data$freedomhouse_1 == 0)
NFree.sub

###
specify_decimal <- function(l,m) format(round(l,m), nsmall=m)

#setEPS()
#postscript("All_Spending_Histogram.eps")
pdf("Policy_Change_Histogram_NotFree.pdf")
#tiff("Policy_Change_Histogram.tiff", height=6, width=6, units="in", res=600)
x <- NFree.sub$Percent_changes_in_GRI
View(x)
kurtosis <- specify_decimal(kurtosis(x), 2)
kurtosis
l.moment <- round(samlmu(x), digits=2)
l.moment
l.kurtosis <- l.moment[4]
h <- hist(x, breaks=50, main="Not Free Countries", xlab="Percent Changes in GRI (Weekly)", ylab = "Frequency", ylim=c(0,5000), xlim=c(-.6,.8), col="gray50", border="gray50")
xfit <- seq(min(x), max(x), length=100)
yfit <- dnorm(xfit,  mean=mean(x), sd=sd(x))
yfit <- yfit*diff(h$mids[1:2]*length(x))
lines(xfit, yfit, col="black", lwd=2, lty=1)
abline(v=0, col="black", lwd=2)
abline(v=median(x), col="black", lwd=2, lty=2)
kurt.lab <- bquote(bold(L-kurtosis == .(specify_decimal(l.kurtosis,2))))
med.lab <- bquote(bold(Median == .(specify_decimal(median(x),2))))
text(.5, 4300, labels=kurt.lab)
text(.5, 3950, labels=med.lab)
dev.off()

#############################################################################
############################# Figure 6a: Free & SI Weekly ####################
#############################################################################


##Free

Free.sub <- subset(data, data$freedomhouse_1 == 1)
Free.sub

###
specify_decimal <- function(l,m) format(round(l,m), nsmall=m)

#setEPS()
#postscript("All_Spending_Histogram.eps")
pdf("Policy_Change_Histogram_Free_SI.pdf")
#tiff("Policy_Change_Histogram.tiff", height=6, width=6, units="in", res=600)
x <- Free.sub$Percent_changes_in_SI
View(x)
kurtosis <- specify_decimal(kurtosis(x), 2)
kurtosis
l.moment <- round(samlmu(x), digits=2)
l.moment
l.kurtosis <- l.moment[4]
h <- hist(x, breaks=50, main="Free Countries", xlab="Percent Changes in SI (Weekly)", ylab = "Frequency", ylim=c(0,8000), xlim=c(-1,1), col="gray50", border="gray50")
xfit <- seq(min(x), max(x), length=100)
yfit <- dnorm(xfit,  mean=mean(x), sd=sd(x))
yfit <- yfit*diff(h$mids[1:2]*length(x))
lines(xfit, yfit, col="black", lwd=2, lty=1)
abline(v=0, col="black", lwd=2)
abline(v=median(x), col="black", lwd=2, lty=2)
kurt.lab <- bquote(bold(L-kurtosis == .(specify_decimal(l.kurtosis,2))))
med.lab <- bquote(bold(Median == .(specify_decimal(median(x),2))))
text(.6, 6850, labels=kurt.lab)
text(.6, 6300, labels=med.lab)
dev.off()

#############################################################################
############################# Figure 6b: Partly Free & SI Weekly #############
#############################################################################


##Partly Free

PFree.sub <- subset(data, data$freedomhouse_1 == .5)
PFree.sub

###
specify_decimal <- function(l,m) format(round(l,m), nsmall=m)

#setEPS()
#postscript("All_Spending_Histogram.eps")
pdf("Policy_Change_Histogram_PartlyFree_SI.pdf")
#tiff("Policy_Change_Histogram.tiff", height=6, width=6, units="in", res=600)
x <- PFree.sub$Percent_changes_in_SI
View(x)
kurtosis <- specify_decimal(kurtosis(x), 2)
kurtosis
l.moment <- round(samlmu(x), digits=2)
l.moment
l.kurtosis <- l.moment[4]
h <- hist(x, breaks=50, main="Partly Free Countries", xlab="Percent Changes in SI (Weekly)", ylab = "Frequency", ylim=c(0,7000), xlim=c(-1,1), col="gray50", border="gray50")
xfit <- seq(min(x), max(x), length=100)
yfit <- dnorm(xfit,  mean=mean(x), sd=sd(x))
yfit <- yfit*diff(h$mids[1:2]*length(x))
lines(xfit, yfit, col="black", lwd=2, lty=1)
abline(v=0, col="black", lwd=2)
abline(v=median(x), col="black", lwd=2, lty=2)
kurt.lab <- bquote(bold(L-kurtosis == .(specify_decimal(l.kurtosis,2))))
med.lab <- bquote(bold(Median == .(specify_decimal(median(x),2))))
text(.6, 6000, labels=kurt.lab)
text(.6, 5500, labels=med.lab)
dev.off()

#############################################################################
############################# Figure 6c: Not Free & SI Weekly ################
#############################################################################


##Not Free

NFree.sub <- subset(data, data$freedomhouse_1 == 0)
NFree.sub

###
specify_decimal <- function(l,m) format(round(l,m), nsmall=m)

#setEPS()
#postscript("All_Spending_Histogram.eps")
pdf("Policy_Change_Histogram_NotFree_SI.pdf")
#tiff("Policy_Change_Histogram.tiff", height=6, width=6, units="in", res=600)
x <- NFree.sub$Percent_changes_in_SI
View(x)
kurtosis <- specify_decimal(kurtosis(x), 2)
kurtosis
l.moment <- round(samlmu(x), digits=2)
l.moment
l.kurtosis <- l.moment[4]
h <- hist(x, breaks=60, main="Not Free Countries", xlab="Percent Changes in SI (Weekly)", ylab = "Frequency", ylim=c(0,6000), xlim=c(-1,1), col="gray50", border="gray50")
xfit <- seq(min(x), max(x), length=100)
yfit <- dnorm(xfit,  mean=mean(x), sd=sd(x))
yfit <- yfit*diff(h$mids[1:2]*length(x))
lines(xfit, yfit, col="black", lwd=2, lty=1)
abline(v=0, col="black", lwd=2)
abline(v=median(x), col="black", lwd=2, lty=2)
kurt.lab <- bquote(bold(L-kurtosis == .(specify_decimal(l.kurtosis,2))))
med.lab <- bquote(bold(Median == .(specify_decimal(median(x),2))))
text(.6, 5150, labels=kurt.lab)
text(.6, 4700, labels=med.lab)
dev.off()



